Numerical stability of a new conformal-traceless 3+1 formulation 

of the Einstein equation 



Pablo Laguna and Deirdre Shoemaker 
Center for Gravitational Physics and Geometry, and 

Center for Gravitational Wave Physics 
Penn State University, University Park, PA 16802 

There is strong evidence indicating that the particular form used to recast the Einstein equation 
as a 3+1 set of evolution equations has a fundamental impact on the stability properties of numerical 
evolutions involving black holes and/or neutron stars. Presently, the longest lived evolutions have 
been obtained using a parametrized hyperbolic system developed by Kidder, Scheel and Teukol- 
sky or a conformal-traceless system introduced by Baumgarte, Shapiro, Shibata and Nakamura. We 
present a new conformal-traceless system. While this new system has some elements in common with 
the Baumgarte-Shapiro-Shibata-Nakamura system, it differs in both the type of conformal trans- 
formations and how the non-linear terms involving the extrinsic curvature are handled. We show 
results from 3D numerical evolutions of a single, non-rotating black hole in which we demonstrate 
that this new system yields a significant improvement in the life-time of the simulations. 
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00 ■ So far, every 3D effort to perform 3+1 numerical simulations of space-times containing black-hole singularities, 
or compact objects such as neutron stars, has encountered the frustration of producing short-lived evolutions. The 
impediment is due to the presence of sever instabilities. There is a common agreement among researchers that the 
origin of these instabilities is not purely numerical. Inconsistent boundary conditions far away from the holes or 
' neutron stars, unsuitable choices of gauge or coordinate conditions, inappropriate algorithms to excise the black hole 
singularity, etc. are all contributing factors. Since approximately the mid 90's, the numerical relativity community 
has been accumulating evidence supporting the view that a crucial aspect determining the life-time of simulations is 
the form of the Einstein equation that is implemented in the numerical codes. 

When viewed as an initial value problem, the Einstein equation can be cast in what is commonly known as a 3+1 
form; that is, the Einstein equation becomes a set of constraint and evolution equations that govern the history of the 
geometry of space-like hypersurfaces, geometrodynamics. There are many ways of writing the Einstein equation in a 
^ | 3+1 form. During the early years of numerical relativity, the most popular form was that introduced by Arnowitt, 
I \ Deser and Misner (ADM) |IJ. It is safe to say that all of the attempts to obtain long-term stable and convergent, 3D 
£~ I ■ simulations of black holes using codes based on the ADM formulation have failed. These efforts have produced short 
lived evolutions lasting < 200M, with M the mass of the black hole. Unfortunately, the reasons behind this failure 
^ i are still not known. Nonetheless, because of all of these failed attempts, the numerical community has practically 
abandoned the ADM system. 

One alternative to ADM has been formulations of the Einstein equation that are explicitly hyperbolic gj . These 
formulations have been introduced in part to facilitate handling both the outer boundary of the computational domain 
and the excision boundary in the vicinity of the black hole singularity. Hyperbolic formulations certainly have the 
advantage that, by knowing the characteristics of each of the evolved fields, one is able to apply the correct boundary 
conditions to those fields that require them. Hyperbolic formulations also enjoy the advantage of the mathematical 
machinery to prove existence, uniqueness and well-posedness. Unfortunately, hyperbolic formulations have so far been 
able to produce evolutions of single black holes lasting ~ 600Af — 1300M 0. 

Another formulation of the Einstein equation that has been relatively successful is the one re-introduced by Baum- 
garte and Shapiro and originally developed by Shibata and Nakamura 0]. This formulation, commonly known as 
the BSSN system, belongs to the class of conformal traceless formulations (see next section for a detailed derivation). 
There are several ingredients that go into the derivation of these equations. Besides the use of conformal transforma- 
tions and a traceless decomposition of the extrinsic curvature, perhaps the most important aspect is the introduction 
of a connection. With this connection, the formulation contains some hyperbolic flavor. Studies have shown, however, 
that 3D, single black hole evolutions are still limited to life-times ~ 500M 

Because of the relative success of the hyperbolic and BSSN systems over the ADM formulation, there seems to 
be a certain degree of consensus that a successful formulation should have some hyperbolic flavor. The formulation 
we introduce in this paper retains this property. The innovative aspect of this system enters in shifting the focus 
from the principal part of the evolutions equations, i.e. those terms containing differentiation, to non-linear algebraic 
terms. Our work is motivated by a recent numerical study of evolutions of single black holes in spherical symmetry 
JjJ using the ADM system. In Ref. 0, long-term stable and convergent evolutions were achieved by "adjusting" the 
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ADM system. The adjusting consisted of eliminating or switching the sign of non-linear terms involving the extrinsic 
curvature. The stability that was obtained was remarkable. Unfortunately, an extension to 3D is not possible because 
the adjustment took advantage of the spherical symmetry in the problem. As we will show, it is possible to achieve a 
similar elimination or switching of sign of non-linear terms. Instead of working with the ADM system, one needs to 
consider a modified form of the BSSN formulation. We will present 3D numerical simulations of single, non-rotating 
black holes in which we demonstrate a remarkable improvement in the stability and duration of the evolutions when 
compared with those based on BSSN. 



II. NEW CONFORMAL TRACELESS EQUATIONS 



In order to gain a better understanding of the elements present in the conformal traceless system here considered, 
we will review the derivation of the BSSN formulation. The traditional derivation of the BSSN system starts with 
the ADM formulation. The ADM primary variables are the spatial metric, g^ } and extrinsic curvature, Kij, of the 
space-like hypersurfaces. The equations that govern the dynamics of g^ and are 

doQij = -2aK i:j (1) 
d Kij = — ViVjQ! + a Rij 

+ a(KK ij -2K ik K k j ), (2) 

where d = dt — £p and K — gijK lJ . Although we concentrate our attention on the vacuum case, the inclusion of 
matter sources is straightforward. Above, V-j and Rij denote covariant differentiation and Ricci tensor associated 
with <7jj, respectively. The lapse function, a, and shift vector, /3 l , encapsulate the diffeomorphism freedom of General 
Relativity and are thus freely specifiable. In terms of the ADM variables, the Einstein constraints read 

R + K 2 - K i: jK lj = (3) 
VjK ij - V l K = . (4) 

As customary in numerical relativity, evolutions are performed without ensuring that the evolved data continue to 
satisfy the constraints. Eqs. (B) and (^) are only used to monitor the accuracy and stability of the simulations. 
The first step in obtaining the BSSN formulation is to abandon g^ and Kij as primary variables and work instead 

with $, gij, K and Aij. The relationships between these and the ADM variables are 

$ = 1 lng 1 / 2 (5) 
9ij = e~ 4<E gij (6) 

p -4* A.. 

ly — C .fly , 



— e A^ , (7) 



where A+j = Kij — gij K/3. The first point to notice is that (||) and (||) imply that the conformal metric g^ has 
unit determinant. Enforcing this condition in numerical evolutions has an important impact on the stability of the 
simulations. Unless noted, indices will be raised and lowered with the conformal metric gij. Given (0)-(@), it is not 
difficult to show that the ADM system (|j])-(p|) of evolution equations take the form 

= -\aK (8) 
6 

9o9ij = -2 a A^ (9) 

d K = -ViV*a + a (iyi" +K 2 /3) (10) 
doAij = er^ {-V l V J a + aR lJ ) TF 

+ a{KA ij -2AiiA l j ). (11) 

Before continuing, there are several points to be made about the derivation and structure of these equations. First, 
in deriving Eq. (|l0|), the Hamiltonian constraint (|^) was used to eliminate the Ricci scalar. In Eq. (|ll]), the superscript 
TF denotes the trace-free part of the tensor between brackets, e.g. T? F = T\j — gijg kl Tki/3 for any tensor Tij. When 
calculating Rfj F , the trace R is again eliminated using the Hamiltonian constraint (||). To simplify notation, we did 
not explicitly transform the terms involving covariant derivatives of the lapse or the Ricci tensor, Rij. These terms 
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should be understood as being given respectively by: 

+ 2s y V'*Vja (12) 

+ 4Vi*V i *-4flyV , *V|$, (13) 

where is the Ricci tensor computed from g^ . One should also keep in mind that g^ and Aij are tensor densities 
of weight —2/3, so their Lie derivatives are 

C p = p k d k A lJ + A lk d 3 p k + A kj d t f3 k - 2 -A V} d k {3 k (14) 
and a similar expression for cjij. The function $ is not a true scalar, see Eq. (|5|). Its Lie derivative is given by 

Cp $ = f3 k d k $> + ~ <S>d k f3 k . (15) 
o 

Next is perhaps the key ingredient of the BSSN system. One introduces a conformal connection 

I ' U ji r, <),,■'. (16) 

where the T % > k are the connection coefficients associated with g^j. The last equality holds because of the condition 
9 = 1. The motivation for introducing these connections is to substitute in the Ricci tensor Rij all the derivatives of 
9^ k ^jk m f avor °f derivatives of F*. In doing so, the Ricci tensor yields a system of evolution equations |i[ |{| with 
"hyperbolic flavor." 

Expanding the number of primary variables in the system to include the conformal connections implies the need of 
an additional evolution equation. From definition (|l6|), one can show that the evolution equation for T l is given by 

+ 2at) k A jk + 12 aAVdj® - -a^K . (17) 

o 

In deriving this equation, the momentum constraint (^) was used to eliminate derivatives of A l i . A final point to be 
made, T l is a vector density of weight 2/3, so 

£ p r = [3 k d k r - f k d k p l + 1 f i d k f3 k . (is) 

In general terms, the BSSN system involves: conformal transformations, a trace-free decomposition of the extrinsic 
curvature, the introduction of a conformal connection and the use of constraints to eliminate the Ricci scalar and 
derivatives of the trace- free extrinsic curvature. 

We are now prepared to derive the new system of evolution equations. Because this system retains the main 
ingredients of the BSSN formulation, its derivation proceeds along a similar path. However, the starting point is 
modified. Instead of we use 

(19) 

(20) 
(21) 
(22) 
(23) 

with n a parameter to be determined later. Notice the main differences. We allow for the trace of the extrinsic 
curvature and the lapse function to be conformally transformed. K and N are densities. Also, we use as primary 
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variable the mixed-index form of the trace-less part of the extrinsic curvature, A % j. If n = 0, one recovers the BSSN 
scalings. With the above transformations, the evolution equations take the following form 

d ® = -Ink (24) 

dofa = -2NAa (25) 



d a K 



„6n$i 

Vjv i 

'2 , 



+ (1 -3ra) NKy3 (26) 
d A l 3 = e 6n * (-V i V j a + aR i j ) TF 

+ (1 -n)N KA l j (27) 

i)„V = & k d ik F + \g ij d jk p k -2A»djN 

+ 2NV) k A jk - -Ng ij djk 
3 

+ 12(1 - n)NA ij d 3 <&. (28) 

Here again, the terms involving a and Rij should be understood as being given by dl2j ) and (|l3|), with the index 
raised with . In addition, the lapse expression ([l2]) should be modified with the appropriate substitution in terms 
of the densitized lapse N. 

By fixing the attention to Eq. ( |27| ) , one can understand the motivation behind the particular form of the conformal 
scalings and the use of A 1 j as a primary variable. If we choose n > 1, the term N KA l j in Eq. ( p7f ) could be 
eliminated or have its sign switched. In the spherically symmetric numerical study in Ref. M, it was demonstrated 
that this algebraic term was responsible for the instabilities of the ADM-based simulations. In Rcf. the changes 
in sign or eliminations of these terms were obtained by substitution of the momentum constraint. This adjustment of 
the ADM system is not possible in non-spherical symmetry. However, it is clear that the conformal scalings we have 
introduced accomplish the same effect. 

The importance of eliminating N KA l j can be understood by the following hand- waving argument. Equation (|27j ) 
has the following structure 

d t A l 3 = CA l 3 

+ OA Terms + Terms without A (29) 

where 

C = nd k k + (l-n)N K . (30) 

At the continuum level, for a single, static black hole, the r.h.s. of Eq. (|2^) should vanish identically. However, 
because of truncation errors, this will not be the case in numerical evolutions. If this cancellation does not take place, 
a solution of the form A 1 j oc e ct is in principle allowed. Thus, C > has the potential of yielding an exponential 
growth. Almost all the numerical evolutions of single black holes have been performed with space-time foliations 
such that K > 0. Since the simulations of physical interested are future directed (i.e. a > 0), then N K > 0. As 
a consequence, the second term in C has the "wrong" sign. Because of the conformal scalings we have introduced, 
values of n > 1 correct this sign problem. The first term in C gives the impression that it also has the wrong sign. 
This is not the case. All of the single black hole solutions of interest have shift vectors for which d k (i k < 0. For more 
generic, multiple black hole situations, it should be possible to construct shift vectors that preserve this property. 
If not, a large enough value of n should be able to yield the appropriate sign of C. Finally, we point out that this 
hand-waving argument provides a reasonable explanation of why the simulations in || with K — const < were 
significantly more stable than those with positive K. 

The use of a mixed-index, traceless extrinsic curvature in 3D studies is not new. One of us used this approach in 
numerical studies of inhomogeneous inflationary space-times |)| . The evolution equations used in || did not include a 
conformal connection T % nor a densitized lapse. The use of densitized variables had a different motivation in ||. With 
those densitized variables, the shift terms in the evolution equations could be written in an almost flux-conservative 
form. However, it was quickly realized that, because during an inflationary epoch the determinant of gij grows 
exponentially, working with positive weight, densitized variables becomes eventually problematic. Our attempts to 
obtain stable 3D, single black hole evolutions with a code based on the evolution system in Ref. |l have failed. Given 
the success of the results in the present work, we believe that the failure of the system in 0] is likely due to the 
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absence of the conformal connection T\ There are indications that in spherical symmetry, using flux-conservative 
equations seems to extend the life of the simulations JhJ . This view is also supported by recent work using a hyperbolic 
formulation [ pT| in which the evolution equations are written in a flux-conservative form. Finally, mixed-index extrinsic 
curvature, not traceless, as a primary variable has been used in ID evolutions of plane-symmetric cosmologies 
and gravitational waves p"3j, supporting the view that this form is more advantageous. 



III. STABILITY EXPERIMENTS 



To test the stability of the new system, we concentrate our attention on a single, non-rotating black hole. We use 
the ingoing-Eddington-Finkelstein form of the black hole solution. The black hole singularity is handled via excision. 
That is, a region inside the horizon, containing the black hole singularity, is removed from the computational domain. 
Our excision has a cubical shape. At its boundary, the values of the primary variables ($, gij, K, A % j, T 1 ) are 
obtained by third order extrapolation of the corresponding evolved values. At the outer boundary, the evolved data 
are blended to the analytic solution. Temporal updating is performed via the iterated Crank-Nicholson method. After 
each time-step, we impose the constraints g = 1 and A 1 » = as follows: 

A* j -> A l ^g^A k k /3 (31) 
9ij -» 9~ 1/3 9v- (32) 

Once the new g^ is obtained, we use definition ( |ilf ) to also reset T 1 . 

We do not write the shift terms in an almost flux-conservative form as done in Ref. jjj. All the discretization of 
spatial operators is done with centered finite differences. The exception is shift terms of the form (3 k dk- Those terms 
are approximated using an up- wind discretization described in Ref. 

Since it is well known that the choice of gauge or coordinate conditions has a fundamental impact on the duration 
of the evolutions, we decided to test the new system on the most challenging case, namely the gauge condition in 
which N and (3 1 are obtained from the exact analytic solution. 

Figure [l] shows the L2 norm of the Hamiltonian constraint and the error in the mass of the black hole from a 3D 
evolution. The location of the outer boundary of the computational domain was set to 9M and the grid-spacing to 
0.25M. We imposed quadrant symmetry, namely reflection symmetry across the x = and y = planes. This choice 
was motivated by the work in Ref. Q, in which long term stability was only achieved in octant symmetry and with 
live gauge conditions. We performed evolutions with higher resolutions in octant symmetry only for the purpose of 
testing second order convergence in our code. An identical, quadrant-symmetric, run using the BSSN system crashed 
at ~ 75 M. Short-lived, BSSN evolutions with analytic lapse and shift were also reported in Ref. ||. As one can see 
from Fig |l|, the new system is able to extend the life-time of the simulation an order of magnitude beyond that of 
BSSN. To test the effects of the symmetry conditions, we repeated a simulation using bitant or equatorial symmetry. 
Figure |^ shows the outcome of the simulation. A comparison of the results in Figs. |l] and || indicate that the duration 
and quality of the data still depend on the symmetries imposed. 

Although we have managed to extend the duration of the simulation, as seen in Fig. ^, there is a late-onset of 
an instability. Preliminary experiments indicate that the instability is likely to be connected to the type of outer 
boundary condition and its placement, as well as the symmetries imposed. Our results are, nonetheless, comparable 
to evolutions in Ref. H based on a parametrized hyperbolic formulation. In that study, evolutions lasting ~ 600M 
were obtained for similar computational domain sizes, but without any symmetries applied, using analytic densitized 
lapse and shift conditions. A direct comparison with the results in || is not straightforward due to the numerical 
technique employed, namely spectral methods. 

Finally, we performed evolutions in octant symmetry for 1 < n < 2. We obtained similar stability to that with n = 1. 
The trend observed in these simulations suggests that for values n > 2 it is very likely to experience a deterioration 
on the stability of the evolutions. This should not be surprising. Looking at the conformal transformations (pl|)-(p3|), 
one sees that large values of n have the effect of increasing the spatial gradients of these primary variables. 



IV. CONCLUSIONS 



This paper introduces a conformal-traceless formulation of the Einstein equation that yields a substantial improve- 
ment in the duration of evolutions. This system has various elements in common with the BSSN system. However, 
it has the following new aspects: (1) modification of the conformal scalings of the extrinsic curvature, (2) use of the 
mixed-index form of the traceless extrinsic curvature as a primary variable, and (3) densitization of the lapse function. 
The reason behind these modifications is to eliminate or switch the sign of non-linear terms involving the extrinsic 
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curvature. Our results support the view that these terms are partially responsible for rendering the evolutions unsta- 
ble. We are currently investigating the degree to which each of the above elements of the formulation here introduced 
impacts the stability of the simulations |3 . We have demonstrated that our 3D simulations of a single black hole 
exhibit a remarkable improvement when compared with those using the BSSN formulation. This is especially apparent 
when one considers we have performed these tests in what may be the most troublesome gauge choice, namely ob- 
tained directly from the analytic black hole solution. Furthermore, we did not impose outgoing boundary conditions. 
These conditions have been shown to improve the duration of simulations in hyperbolic and BSSN systems. We are 
currently extending our study to investigate whether the stability of this new system is preserved and improved using 
live gauge conditions and in dynamical situations such as distorted and binary black holes. 



V. ACKNOWLEDGMENTS 



This work was supported by NSF grants PHY-9800973 and PHY-0114375. We want to thank A. Ashtekar for 
comments and helpful discussions. We also want to thank the rest the MAYA Project development team: D. Fiske, 
B. Kelly, E. Schnetter and K. Smith. 



[1] R. Arnowitt, S. Deser and C.W. Misner, in Gravitation: An Introduction to Current Research, edited by L. Witten (John 
Wiley, New York, 1962), pp. 227-265. 

[2] For reviews on hyperbolic methods for Einstein's equations: "The Cauchy problem for the Einstein equations" , H. Friedrich 
and A. Rendall, in Einstein's Field Equations and their Physical Interpretation, ed. by B. G. Schmidt, Springer, Berlin, 
2000. Published in Led. Notes Phys. 540, 127 (2000). e-Pri nt Archive: gr-qc/000207. Al so: "Hyperbolic methods for 
Einstein's equations", O. Reula, Living Reviews in Relativity (http://www.livingreviews.org) 

[3] Lawrence E. Kidder, Mark A. Scheel, Saul A. Teukolsky, Phys. Rev. D, 64 064017 (2001) 

[4] T. Baumgarte and S. Shapiro, Phys. Rev. D, 59 024007 (1999). 

[5] M. Shibata and T. Nakamura, Phys. Rev. D, 52, 5428 (1995). 

[6] M. Alcubierre and B. Bruegmann, Phys. Rev. D, 63 104006 (2001). 

[7] B. Kelly, P. Laguna, K. Lockitch, J. Pullin, E. Schnetter, D. Shoemaker and M. Tiglio Phys. Rev. D, 64 084013 (2001). 

[8] A. P. Gentle, D.E. Holz, A. Kheyfets, P. Laguna, W.A. Miller, D.M. Shoemaker Phys. Rev. D, 63 064024 (2001) 

[9] H. Kurki-Suonio, P. Laguna and R.A. Matzner Phys. Rev. D, 48, 3611 (1993) 
[10] L. Lehner, private communication (2001). 
[11] D. Neilsen, in preparation (2002). 
[12] J. Centrella and J. Wilson, Ap.J., 273, 428 (1983 ). 
[13] J.M. Bardeen and L.T. Buchman, |gr-qc/0111o"85 . 
[14] A. Ashtekar, P. Laguna and D. Shoemaker, in preparation (2002). 



7 



3.3 
3.25 

CV2 

o 3.2 



y 3.15 
3.1 



2 


2 



CO 

I 

o 



S -4 



-6 



i i i i i i i i i i i i i i i i i i r 



_l i i i l i i i l i i i l i i i l i i L 

200 400 600 800 

t/M 



FIG. 1: L2 norm as a function of time of the Hamiltonian constraint (HC) and the error in the mass of the black hole (M — 1) 
from a 3D evolution of a single, non-rotating black hole in ingoing-Eddington-Finkelstein coordinates. The outer boundary of 
the computational domain was located at 9M with M the mass of the black hole in the initial data. The grid-point spacing in 
this run was 0.25M. The run was performed imposing reflection symmetries across the x = and y = planes. 




FIG. 2: Same as in Fig. y but imposing reflection symmetry only across the x — plane. 



